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1.0 Summary 


This report contains the development of a recursive, non-singular method for computing 
the first and second partials of the gravitation potential, with respect to position, using 
both unnormalized and normalized harmonic coefficients. When unnormalized coeffi- 
cients are used, every attempt was been made to build a "fast” algorithm. When normal- 
ized coefficients were used, the algorithm developed uses a more stable, albeit more 
complex, recursive algorithm for the derived Legendre functions. Even so, the normal- 
ized algorithm is still quite efficient. The normalized algorithm should be quite stable 
and portable for model sizes exceeding 180x180 in degree and order. Efficiency in com- 
putation was gained by precomputing everything that was not a function of the state and 
by using singly dimensioned arrays wherever possible as well as arrays of pointers to 
arrays. 

A complete derivation of the gravity gradient torque resulting from a full (nxn) gravity 
model is given since it uses the second partial of the potential developed earlier. 

A complete derivation of the geomagnetic field vector was Included since the computa- 
tion of the magnetic field is so similar to that of the gravitational field. 

Ada code for all of the algorithms is included. 

Test cases compare the algorithms to each other and to previously published data. 
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2.0 Introduction 

This report is basically a rewrite of Ref |2J, with some useful additions. First of all. by 
examining the derived Legendre functions that are used to compute the gravitational 
acceleration it is noted that some of them sue not functions of the state and hence may 
be computed only once. This fact is used to speed up the computation of gravity and its 
partials. 

Secondly, a derivation using normalized gravity coefficients and a superior recursion for- 
mula for the derived Legendre functions is presented. As the size of gravity models 
increases, an algorithm using normalized coefficients becomes more attractive since the 
unnormalizing process requires the computation of terms on the order of 2n! Even for 
models of size 50x50 this would be a number so large ( - 10 158 ) that some computers 
might not be able to compute it. 

Algorithms are developed, using both normalized and unnormalized gravity coefficients, 
that compute the first and second derivatives of the potential function. This yields the 
gravitational acceleration, and the partial derivative of the gravitational acceleration with 
respect to the position vector. The partial derivative of the gravitational acceleration is 
needed in the computation of the state transition matrix for both estimation and optimi- 
zation. In addition, the partial derivative matrix can be applied to the problem of com- 
puting general gravity gradient torque. 

Next, a general gravity gradient torque derivation is presented that uses the second par- 
tial of the potential developed in the previous section. 

Since the geomagnetic field is defined in terms of Legendre functions, a derivation of the 
geomagnetic field is Included which is very similar in form to the gravity derivation. 

And finally, Ada code as Implemented in the Ada Simulation Development System 
(ASDS), [9J. is given for the various algorithms in addition to test cases that verify the 
validity of the derivations. 
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3.0 The Gravitational Potential Function 

The gravitational potential function is normally written 


V= "7"£ P r*.m( e H C n.mCOS(mX)+S n>m sin(mX)) (3-1) 

n = 2 m = 0 

where ji is the gravitational constant, a e is the equatorial radius, r is the magnitude of 
the position vector, X = (x v x%, x 3 ) , and 






(3-2) 


are the associated Legendre functions and P n are the Legendre polynomials. Also, we 
have the sine of the latitude 


e = x,/r 


and the longitude is computed from 


(3-3) 


X, 

tanX = — 

For notatlonal convention, we define a potential function U to be 

17* -V 

and write eq(2-l) as 


(3-4) 


U= 7 + S Z7 ( T C) p n.m(e)(C n>m cos(mX) + S nifn sin(mX)) 


n-2m-0 


(3-5) 


Given the equation for e above, becomes 






3e m r" 1 


where 


(3-6) 


p 2 *X? + j| 

P>— " 

" de m 

The are known as derived Legendre functions. 


(3-7) 


Note that 
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( 3 - 8 ) 


dp 

*1 n 

d 2 f 

LP P J 

d\ 

~*2 *1 n 


.P P J 


Also define. 


C m sp m cosmX 

S m *p m sinmX ( 3 - 9 ) 

B n ,m mC n,m C m+ S n,m S m 

The m and the m are the unnormalized cosine and sine gravity coefficients that 
result from the mass distribution of the planet. When these coefficients are published, 
they are published In normalized form. The relationship between the normalized and 
unnormalized form Is given by 


C am* W ( n ’ m ) C n.m 

B n,m~ N(ru m) S n%m 

where 


( 3 - 10 ) 


N(n, m) = 


(n-m)!(2n+l)(2-6 om ) 
. (n+m)! 


1/2 


( 3 - 11 ) 


Where 5„ m is 1.0 If m = 0, and is zero otherwise. 

The derivation will proceed using unnormalized notation because the derivation is some- 
what simpler. In a later section, a derivation using normalized coefficients will be devel- 
oped. The normalized form does not require the computation of terms on the order of 
(n+m)! This may be desirable on some computers where very large or small numbers may 
cause a problem. For now, the potential can be written 



ii a "P 1 ” 

ns 2ms 0 


( 3 - 12 ) 


This form Is especially useful since C m , and S m can be calculated recursively and 
the singularity at the pole (p = 0) can be avoided. 


The unnormalized derived Legendre functions may be calculated recursively a number of 
different ways. In [1] and 12] the were computed from 


C = C-2+ (2n-l)Ci 1 .(m^D 
p£ = P„ = ( (2n- l)eP„_, - (n- 1)P „_ 2 )/n 
Pg= 1 f^ = o 

P? = e P] = 1 
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In [3], seven recursion algorithms were compared numerically for stability. Unfortu- 
nately, eq(3-13) was not among those studied. Of the seven algorithms studied, two 
were clearly superior. The simpler of these is 

P^= ((2n-l)eP^_i- (n+m-l)P^. 2 )/(n-m), (m<n) 0-14) 

Note that when m=0, eq(3-14) reduces to eq(3-13) for m=0. Experiments similar to those 
carried out in [31 were conducted by the author using the normalized error between a 
single precision computation of the Pjj* and a double precision computation of the 
using both eq(3-13) and eq(3-14). In every case eq(3-14) had lower error. The worst error 
in all cases occurred for e = 0.2. 

Although eq(3-14) is highly stable, it cannot generate the diagonal elements Pj|. Realizing 
that Pj |_ 2 = 0, (all PJJ 1 beyond the diagonal are zero), eq(3-l3) can be used to compute 

K = (2n-l)P^l} (3-15) 

Starting with either eq(3-13) or eq(3-14) it is rather easy to show that the inner diagonal 
terms, PJJ " 1 , can be computed from 


P" -1 = eP" 

* n cl n 

Also, note that 


(3-16) 


C m = P m cosmA, = p m " 1 + 1 cos(m- 1+ 1)X = CjC m _i - SjS m _, 
S m = p m sinmX = p m - 1 + 1 sin (m- 1 + 1 ) X = S, C m _ , + C, S m _ , 


JC, 

C, = pcosX = p— = x, 


(3-17) 


S, = pslnX = p— = Xj 

It is helpful during coding to note thatC m e C m /x ,n andS m = S /r m are also recursive, since 




Cj C m _ i - Sj S m _ i 


Si C m _ i + Ci S m _ i 



(3-18) 
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4,P The First Partial 

The gravitational acceleration vector, g, is calculated as the first partial derivative of U 
with respect to the planet-fixed position vector, & From eq(3-12), we have 


where 


_dU_dUdr dU& dU dB n ,m 
9 3y drdX+dthX dB. m dX 


(4-1) 


31/ P T'y' P , a e, " ( n + nt+ 1) 

a? = J2 Z^j 1 r ' r™ n n - m 


(4-2) 


317 

Tt 




n, m 

r^~ 


(4-3) 


at/ 


as 


n, m 


ii a n P m 

lX-r<T> fi 


(4-4) 


We note that 


3Bn,m _ _ 3c m ap ac^x as m ap as m ax 
ay Cn>m( 3p ay + 5X a^ &n - m( 5p ay 5X ay 


or 


(4-5) 


^ n,m = mC^ m (p m_ 1 cosmX^-p m sinrnx|^) +rnS n>m (p m 1 slnmX^+p m cosmX^) 

dp ax 

Now, using the definitions of and 5 — given earlier, 


dB n m 

n>m = mo m_1 C 

3v "t '“'n, m 


ay 


cosmX 


slnmX 

~*2 

\ 

P 

*2 

. 0 . 

P 

x i 

. 0 . 

> 


+ mp m-, S^ m 


slnmX 


cosmX 

"*2 

\ 

P 


+ 

P 

*1 



_ 0 _ 


_ 0 _ 

/ 


Xj 

where column matrices are used in the Interest of saving space. Since cosX = — 
and sinX = , and using the definition of C m and S m . it Is rather easy to show that 


Also, 


3B, 


n, m 


ay 





B n, m- 1 

- 

— m ( C r, m S m _, -S^ m C m _j) 

s m 

~^n, m- 1 


0 


0 


e mb 


(4-6) 


— *• n 

1. Note: For notational simplicity, * X X throughout 

n=2meO 
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(4-7) 


and 


3^ = X 
dX T 


3e 

dX~ 



where 

b t = ( 001 ) 

Combining these partlals and substituting Into eq(4-l), we get 


Defining 


H f, , VV 1 i\ 




a n pm 

(_? ) n 
r ^>-1 




nS'm- 1 + ^n, m®m- l) 


0 


n, mPm- 1 ) 


n pm 

^n = m ^TT ( C n, m C m- 1 + S n, m S m- 1 ) 
m* 1 

n pm 

m p^rj ( C n. m S m- 1 ~ S n. m C m- 1 ) 

m* 1 

p™ 


r „=C ni0 (n + l)P^ + 2 ( 1 + n+ m) - j-j (C am C m + m S m ) 


m* 1 

n pm+ 1 


Wn-C n , 0 Pi+ 2 -V( C n.mC m +S„ >m S m ) 


m* 1 


(4-8) 


(4-9) 


(4-10) 


(4-11) 
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we can write 


_ a n 

Hs X <f> Hn 


n = 2 


_ a n 

'■ 1 + S<T> r » 


n » 2 


_ a " 

j- i < T *> j. 


n = 2 


_ a n 
K= £ <-^> K„ 


n = 2 


A = r+ 




*.* 


( 4 - 12 ) 


9 






( 4 - 13 ) 


Note that the final result for the first partial derivative of the potential Is a rather simple, 
compact, vector equation. 
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5.0 The Second Partial 


3 %J 

Next, we calculate starting with eq(4-l), as It leads directly to a compact symmetric 
notation. ^ 

We have from eq(4-l) 


dU _ dUdr dUde dU dB n,m 

dX ~ drdX + de dX dB ^ m dX 


Thus, 


d\f _ dhjdr dr T d 2 Ude de T d 2 U dB^ m dB l m 
3X 2 ~ d r 2 dXdX de 2 dXdX + dB^J* ** 

d 2 U dr de T de dr T d 2 U p B n , m dr r dr dB^ m ) 
dedr ( dXdX dXdX } drdB a m [dX dX dXdX ) 

d 2 U p B n, mte T de dBl m \ dUd^ dud^t dU d 2 B^ m 

dedB^dX dX + dXdX Jdr d x <l+ ^dX 2+ d B n. m dX! 2 


(5-1) 


The second and cross partlals appearing in eq(5- 1) are 


d\f_ 2p 
dr 2 r 5 + 



"(m+n+l) (m+n+2) 


K B *.m 


d 2 U 

de 2 



dhj 


= 0 


dh; 

dedr 




dedB 




n pm 
* n 


+ 1 


n, m 


4mdS-i r 


d^u _ V'V 1 l 1 / a e N "( m+n+ 1) Din 

drdB^- LL^r' r* 

dr dr T _XX T 
dXdX~ r 2 


(5-2) 
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3e de T 
d&X 


= ^ (Xq. t + aX T ) + ^Q-a 1 

r 6 r 4 r 2 

3\_ = 1( 

3S 2 H r 2 J 


^ • = ^X2f T - A (Xa T + &X T ) - -?/ 


32* r 5 


x 3 

7 


3% 

3* 2 


n, m 


= m(m- 1) 


B n, m-2 -^n,m-2 ® 

“A^ m-2 “®n,m-2 ® 
0 0 


where I Is a 3x3 identity matrix, and 


B n, m-2 E ^"n, m^m-2 + ®n, m^m-2 
•^n, m-2 E ^-'n, m®m-2 — ^n, m^m-2 


We also note the special combinations 




and 


<* B n, m3r T 3r <* B n, m_ 

dX dxd2@x 


m- 1 


-A 


A m- 1 

o 


r +m r [®n,m-l "^n, m- 1 ®] 


= ™(bX T +Xb T ) 


where 


b = 


J n t m- 1 


-A, 


0 


and 


(5-3) 


(5-4) 


(5-5) 


(5-6) 


(5-7) 


3B„, m deJ 3 
dX dX d&X 


j (W T + ffJ? T ) - ^ (bX T + Xb T ) ) 


Putting all these Into eq(5-l) leads to 


(5-8) 
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* 2 U _ H f0i (m+n+2) nrn 2 & t 

m*--? { 2 + LL ( t } (m+n+1) — ^ — OW-^- 

■^XX^ C +2 ^[^ T ^(^ff T +fi^ T )+fiff T ] 

+ ?SI <T» n ^ m [<^ T+ ^ T ) - 7 (MW)] 

- 7 II < T> " * < ( < *» T > ■ 

-p (‘ * II <7*> f-] 

*7ll ( 7> Cy[3^B'-|ui« r *rf)-5|j 


P? 


+ ^XX<t> 


B, 

-A 


'n,m-2 ^n,m- 2 ®l 


-A„ 


L rt» m-2 
0 


*n, m-2 ® 
0 0 


(5-9) 
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If we now define 


Le 2+ X ( T ) X (n+m+ (n+m+2)P ^%r 


n = 2 m = 0 


** rv n n prn+ 2 

s-^ b - 

n «= 2 m = 0 


n n 


X ( T c ) 

n= 2 m= 2 

- a n n 

0^(7) 


n-2 m = 2 

a 


n n 


p *E<f ) X'^’* 1(m+n+I) 


( ^n, m^m-2 + &n, m^m- 2 ) 

y/n-2 

^n,m^m-2 ~ ^n.m^m-2) 
r m "*2 

B 


n, m 


n = 2 m* 0 


n n 


0* X (_ F ) X^ +lm 

n = 2 m= ] 

- a n ” 

2 p : 


+ (^n, m^m - I + ]) 

^-1 

( ^n, m®m- 3 ~ ®n, m^m- 1 ) 


+ 1 

m 


n = 2 m= 1 


a. 


n n 


s* X ( T ) X (m+n+1)P " m 

n* 2 m= 1 

T “"X ( T )n S < m+n+1 > p T m 

n = 2 m = ] 

then with these definitions 


r"*- 1 

( ^n. m^m- 1 ~ ^n. m^m- 1 ) 

r'"" 1 

( m®m- 1 “ &n, m^m- J ) 


i-l 


we can write eq(5- 1) as 


ff = 


9 

R 

0_ 


r- 


s 

T 

0 


(5-10) 


(5-11) 
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(5-12) 


0 - + 7 <*» T+ e^ T ) + ®e T _ 

- ^[7 (2fe T + ft* T ) - 2p?2QC r J + ^(ffff T + fi© 7 ) - ^ (©2f T + X© t )J 

- $ . j* r[,- ^ 1 w> . a,] 


Collecting like terms, we get 


11 "-a 0 

+ 1-Q-JV0 

r |_o 00 




£(t + i« 2 f*apS+r*ai 4 )(S£) 
■7 (M t + p+ w psIi5?T] 


r r 3 r 


* £ { (BS’+BS’} - 5 (B^+XB 7 ) - ( 3gr * 30fT ) + 


JV-QO | 

-a -no | 
_o 0 o| J 


Recalling that 


and defining 


A = (r+eH) 



r 


we have finally 


F*L+e(Me + 2(P+H)) +A 
G = -(Afe + P+H) 
file Eff + T 





(5-13) 


(5-14) 


+ 



-1 

0 


_*r_ 


n- a -n e 1 

-Q -(JV+A) R i 
O R -A J 


(5-15) 


Note that the final result for the second partial derivative of the potential is a rather sim- 
ple, compact, symmetric matrix equation. 
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Big. Cpipptttatipnftl Co n s i derations 

If the derived Legendre functions are arranged in a table such as Table 1 below, it is 
rather easy to see that to the right of the diagonal all terms are zero. Along the diagonal 
only pure numbers appear, and immediately to the left of the diagonal, the diagonal term 
appears multiplied by e . This means that a number of the derived Legendre functions 
that are needed in the computation of gravity can be computed once only and stored, 
and do not need to be computed using the recursion relationships. 

This can save a great deal of time. If in addition, all coefficients in the recursion relation- 
ships are computed and stored as functions of n, another slight savings can be gained. It 
was found that these two taken together save about 15 - 20% of the time normally taken, 
depending on the size of the coefficient array used. 

In Ref [1] and [21, the gravity coefficients were placed in a single array in an attempt to 
avoid the time it takes to manage a two dimensional array on a computer. In the new 
code given in Appendix A, it was found that by having an array of pointers to arrays the 
code is more clear and just as fast. 



m=0 

i 

2 

3 

4 

5 

6 

7 

n = 0 

oh 

II 

t-i 

©SL 

II 

o 

o 

II 

etc. 





i 

ii 

m 

p} = 1 

m 

A 

II 

o 

etc. 




2 

r 2 

= 3c 

CO 

II 

V 

I& 

II 

o 

P% = 0 

etc. 



3 



ff= 15e 

^=15 



etc. 


a 



* 

= 105e 

= 105 

EB 

O 

II 

V 

etc. 

5 



*8 


= 94 5 e 

Pf = 945 

P® = 0 

cn^j 

II 

O 


Table 1. Table of Derived Legendre Functions 


In the code given in Appendix A, the “fast* gravity model takes advantage of all of the 
things mentioned in this section. The normalized model is also "fast" in the same sense. 

Since no difference in the computation of the gravitational acceleration could be detected 
out through a 30x30 model using both eq(3- 13) or eq(3- 14) to generate the PJJ 1 . and since 
eq(3-14) requires more multiplications and additions than eq(3-13), eq(3-13) was used 
for the “fast" model. Since normalization makes more sense as model size, n, increases, 
it was decided to use eq(3-14) in constructing the normalized algorithm. (Section 7). 
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7.Q first and Second Partial* Using Normalized Coefficients 


Following (3) and [4], define the normalized derived Legendre functions, 

P?*N(am)P? (7-1) 

where N(n,m) Is given by eq(3-l 1). Note that this definition means that 


Note also that If m = 0, 


p:c n . m =cc n , m 

^n^n, m = m 


and if m > 1 


1 

N(a 0) = (2n+ l) 5 


(7-2) 


(7-3) 


N(am) = ( 


(n-m)! (2n+ 1)2 5 


(n+ m)! 


) 


(7-4) 


The recursion relationships given earlier for the derived Legendre functions are repro- 
duced below for the sake of convenience 


K=K -2 + (2n- l),or 

K = ((2n-l)eP^_ 1 - (n+m-l)P^_ 2 )/(n-m), (m<n) 

P^ = P n = ((2n- l)eP n _j- (n- l)P n _ 2 )/n 

P^= (2n-l)P^Z} (7-5) 

P" _1 = eP" 

‘n cr n 

pg = 1 = o 

P? = e P} = 1 

Either of the formulae for computing the PJJ 1 may be used. The first is faster, and the sec- 
ond Is more stable numerically. The difference In the resulting algorithm is slight, as 
shall be seen. Taking the definition given In eq(7-l) and applying It to the first of eq(7-5) 
gives 


* - N(nm)p z - (2n - 1)N,n * ™ 


or. 


Pn ~ N(n-2, m y n - 2 ( } 


N(am) 

■v j 


N(n- 1, m- 1) 


(7-7) 


Similarly, 
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^=N(n,m)P^=t ^-}l N(am) P? (n + m- 1) _ N(g m) , 

n (n-m) N(n- 1, m) n_1 (n-m) N(n-2,m) n ~ 2 m< " 

(n-1) N(n,0) 


pP_ N( nO) pP-r^ 2n_ N ( n '°) p> ( n ~ 

i--N(n, 0)^-e — - N(n _ 1>0) ^-i ~ 


K=N(nn)PZ= (2n-l) 


N(ru n) 


N(n- 1, n- 1) 


N(n- 2,0) 

l 

p"- 1 _ ,2n+ 1 2-n-j 

r n - 1 ~ ( — > *n-l 


pP , 

* n-5 


C’ = iV(an-l)P;-^E - N ^ - - n _' 1) p; = E(2n-l) 


2 n 
N(ru n- 1) 
N(n- 1, n- 1) 


(7-8) 


n£2 


P 71 ’ 1 

" n- 1 


N(a n) 

F? = 1 Pj = 0 

P? = 73e P} = 73 

These define the recursion relationships for the p£\ at least symbolically. Computation- 
ally, these relationships can be simplified. Defining 


5(n, m) = 
Tj(n,m) 


N(ru m) 


N(n- 2, m) 
(2 n- l)N(n, m) 


_ / (n-m) (n-m - 1) (2n+ 1)^2 
m2l ~ ( (n+ m) (n+ m- 1) (2n-3) } 


N(n- 1, m- 1) 


m2 


/ 

1 7 


2 (2n+ 1) (2n- 1) 


n+ m) (n+ m- 1) (2-5, 


\ 

0, m- 1 ' / 


and 


(7-9) 


5(n, m) 


(2n- 1) N(rum) 
(n-m) N(n- 1, m) 


= (2n- 1) (2n+ 1) 2 

1 (n- m) ( n+ m) ’ 


il (Tim) 


(n+ m- 1) 

N(ru m) 

(n-m) 

N(n- 2, m) 


(n+ m- 1) (2n+ 1) (n-m- 1) 2 
” *■ (n+ m) (n- m) (2n-3) ^ 


the recursion relationship for P% becomes either 

PJT = £ (n, m) P?_ 2 + tj (n, m) P^J, 1 


(7-10) 


(7-11) 


or 


P^ = ^(n,m)eP^_i-Ti(n.m)P^. 2 (7-12) 

Note that £(n, m) and (n, m) or $ (rum) and i) (a m) are constant functions of n and 
m and need be computed only once. There is no need to use both eq(7-ll) and eq(7-12), 
either one will do. The code for the normalized model contained in Appendix A is based 
on eq(7-12). Another model was developed based on eq(7-ll), but that code is not 
included in Appendix A. The test cases in Appendix B labeled “norm_I" came from the 
model using eq(7-l 1), and test cases labeled “nonn_IT came from the model using eq(7- 
12). The only difference in the code is the use of £ (n, m) and q (n, m) in lieu of % (rum) 
and q (n, m) . and eq(7-l 1) in lieu of eg( 7-12). In either case, everything that follows will 
be exactly the same. Since roughly more multiplications are required if eq(7-12) is 
used, it was anticipated that the “norm_I" model would run faster than the “norm_IT 
model. The difference in run time turned out to be less than the noise in the timer. This 


16 


NASA CR- 188243 


January 1993 


111 


means that “normjl" is the way to go, since It is based on a more stable recursion for- 
mula. 

Going back to eq(7-8) and defining 

(2n- l)lV(aO) V(2n+ 1) (2n- 1) 

' nN(n- 1,0) n 

i (7-13) 

B(n) r (a_ 1 > JV < n -°> _ < n “ l) ( 2n+ V 
’ “ nN(n-2, 0) n ( 2n-3' 

the equations for the normalized derived Legendre polynomials become 

= ea(n)P£_j-p(n)31_2 (7-14) 

Note that a(n) and |3(n) need be computed only once and stored since they are only 
functions of n and are not functions of the state. 

Going back to eq(7-8) and defining the inner diagonal term, 8 (n) , as 

8(n)-(2n-l) 5 ^^liL_P^:} = (2n+l) 2 P£:l (7-15) 

The Inner diagonal, PJJ" 1 can be computed as 


PJJ -1 = e8(n) (7-16) 

Note that 8 (n) need be computed only once and stored since it Is only a function of n 
and Is a not function of the state. The use of 8 (n) will speed up the computation since 
only a single multiply Is required to build the inner diagonal term. 

Looking now at eqs(3-29) given earlier for J n , K n , T n , and H n , note that multiplies 
and S^ m everywhere except In H n . Therefore, Pjj\ C am and S_ m may be replaced 
by P%, C n>m and everywhere except In H n where the terms P^ + ' C IKm and 

1 m appear. Multiplying and dividing by N(n,m) and N(n,m+ 1) leads to 


r/n+ 1 p 

r n ^n, m 

Defining 


N(qm) N(n, m+ 1 ) pm^i c 
N(rum) N(rt,m+ 1 ) " "> 


N(n,TTl) 

N(i% m+ 1) n n>m 


note that 


C (am) 


N(n, m) 
N(n, m+ 1) 


'(n-m) (2-8 0 J (n+m+ 1) 


C(n,0) = ( 


n(n+ 1) 
2 


2 


1 

C(n,m) ) mll = ((n-m) (n+m+ l)) 2 


(7-17) 


(7-18) 


(7-19) 
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These last two may also be computed and stored and then used to compute H n . This will 
allow the computation of gravitational acceleration without the necessity of unnormaliz- 
ing the gravity coefficients. This makes the algorithm more portable, since for laige grav- 
ity models, some computers can’t handle the large numbers involved in the 
unnormalizing process. 

An examination of L, M, N, fl, P, Q. R, S, and Twill show that the only other term needed 
in order to compute the second partlals using normalized coefficients is in the term, M, 
and Involves P^ +2 C^ m and PjJ' +2 S^ m . Consequently, form 

pm + 2 r _ JV(qm) N(n,m+2) N(qm) ^ 2 - 

" "•"* AT(am) JV(am+2) " "• m N(n,m+ 2) " ' 

As before, define 


T (a m) 

Note that 


1 

JV(am) _ ^(n-m) (n-m- 1) (2-5 0m ) (n+m+ 1) (n+m+2)y 


N( a m+2) 


( 2 - 8 , 


’o, m+2 


,) 


r 


(7-21) 


T(aO)M n(n - 1)( '!; 1)<n+2) )^ 

1 

T(a m) ) mal = ( (n- m) (n-m- 1) (n+m+ 1) (n+m+2)) 2 


(7-22) 


Again, these may be computed and stored and used to compute M. This will allow the 
complete computation of the first and second partlals of the potential using normalized 
coefficients. 


One comment, the simplification In Table 1 that allowed the Inner diagonal to be the 
main diagonal multiplied by e is no longer valid. In the normalized case, the inner diago- 
nal term is computed by multiplying 5 (n) by e. 

All the normalized equations have been coded in ASDS [9] and verified against previously 
published (2) test cases, and found to agree exactly. The code for the normalized gravity 
model using the recursion relation from 131 is given in Appendix A. and the test case data 
is given in Appendix B. 
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8.0 Gravity Gradient Torque 

The gravity gradient torque on a spacecraft is derived twice here. The first derivation 
uses a point mass gravity model and assumes that l/|r+ 8?) 3 is approximated by the first 
term of the binomial expansion. 

The second derivation uses a full n x m gravity model and assumes that gravity varies 
linearly about the center of mass, i.e., that gravity in the vicinity of the center of mass is 
given by g = 0 Cfl +|?8r . where both g cg and |f are computed using a general gravity 
model subroutine such as the models discussed in previous sections, and given in 
Appendix A. This derivation shows that the eigenvectors in [7] are not needed, as was 
pointed out in [8]. 

It is then shown that the point mass derivation and the general derivation give identical 
results when the harmonic coefficients of the full gravity model are set to zero. 

It is shown in Appendix B that when only J2 (-Cjo) is used, the general formulation gives 
the same torque as that given by Roithmayr’s (5) model. 

It is anticipated that the use of the full potential model in the calculation of the gravity 
gradient torques will lead to more accurate attitude simulations. 

8.1 Point Mass Gravity Model 

In Figure 1, the vector p describes the position of a particle of mass, dm, in the body axis 
system. 



Figure 1 Position of a Particle of Mass in the Body Axis System 


Assume the matrix, B, relates the body axis system to the system in which r and f 1 are 
defined and g is computed. Then, 


fj = r+Bp 

The gravitational force on dm is 


( 8 - 1 ) 
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( 8 - 2 ) 


dF * -n-idm 

Rotating this force back into the body system yields 


dj = B T dF = -pB T ^dm 


p(B T r+p) 






dm 


We first examine 


r*= 3/2 

= [(f T +p T B T )(f+Bp)] 3/2 


= [ r 2 + p 2 + 2p T B T r] 3/2 

We note that the term B T f appears In both the numerator and denominator 
expression for dj and rewrite It as 


B T r = r B T t 

where r is the unit vector along f . 

Now define the unit vector b In the body axis system 


fe = B T f 

and write 


dj=~ 


p(rb+p) dm 

3/2 


(t 3 + p 2 + 2rp T b) 

The moment about the center of mass due to dj is 


dx = p x df = — 


prp x bdm 


(r 2 + p 2 + 2rp T b) 3/2 


If we factor out r 2 and Ignore p 2 compared to one, we get 

prp x bdm 


dx = — 


, 3/2 


Now, using the binomial theorem and again Ignoring p 2 terms, we get 

dx = ~px£>^l-3^jdm 
Integrating over mass, we get 

x = Jdx = ~Jpdmxb+3^J(pxb) (p’fydm 


(8-3) 

(8-4) 
of the 

(8-5) 

( 8 - 6 ) 

(8-7) 

( 8 - 8 ) 

(8-9) 

( 8 - 10 ) 

( 8 - 11 ) 
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( 8 - 12 ) 


but jpdm = o since p Is measured from the center of mass: hence 

x = ^J(pxfc) ((Fb) dm 

where 


P T f> = Px b I + Py b 2 + P* b 3 

and 


therefore. 


(pxfe) (p T b) 


thus 


t J k 


Py b 3-Px b 2 

Px Py P« 

n 

P* b l - Px b 3 

b\ b 2 b 3 


Px b 2-Py b l_ 


P*Py b l b 3 Py^hj PyP* b 3 
PxP« b ? PyPxhlhS P* b l b 3 

P* b A PxPy b l PxPzhshs 


-PxP* b l b 2 -PyP, b 2 "Px^hs 
“Px b l b 3 ~P*Py b 2 b 3 "PxPz b 3 
~P*Py b ? Py b l b 2 -PyPz b l b 3 


T 


3 ^ ^2 b 3 (I**" fyy) 

"jj b l ^3 (fxx - f**) 
b l b 2 ^fyy ” ficx) 


+ b, bg/jq, -b, bjfj,, 

-b 2 b 3 I xy + b, bj/yx 
+ b 2 ba^xi -b l b 3*yx 


+/ tfz (b2-b2) 
+ fxz( b ?- b 3> 

fxy ( b 2 — h? > 


where we define the moments and products of Inertia to be 


*xx= J(Py + P*)dni 
J vy = J(Px + P^ dm 
I„ = J(P* + Py) dm 
I xy = J(PxPy> dm 
f xz=J<PxP,> dm 
*«* = J<PyP,) dm 


8.2 General Gravity Model 

The force of gravity at the particle, dm. Is now assumed to be computed from 


9 - 


a+^8r 
* c a dr 


(8-13) 


(8-14) 


(8-15) 


(8-16) 


(8-17) 


(8-18) 
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Referring back to Figure 1, we recall that B rotates p back Into the system where g Is 
computed, t.e.. that 


8r = Bp 

and therefore the force on the particle Is given by 

dR = (g CJ7 + |f Bp) dm 

Rotating this force back Into the body system yields 

df = B T dF = (B T g cg + B^Bp) dm 
The moment about the center of mass due to dj is 

dt = p X df = (p X B T g cg + p x B^Bp) dm 
Integrating over mass, we get 

t = Jdx = JpdmxB T g Cfl + JpxB^Bpdm 

The Integral Jpdm = 0 since p Is measured from the center of mass. Hence, 

% = JpxB^Bpdm 

Next, define 

G=B 7 |?Bf 

dr 

then 

x = Jpx Gpdm 


( 8 - 19 ) 


( 8 - 20 ) 


( 8 - 21 ) 


( 8 - 22 ) 


( 8 - 23 ) 


( 8 - 24 ) 


( 8 - 25 ) 


( 8 - 26 ) 


t It must be pointed out that Jjp is computed in an equatorial or planet-fixed system and must be rotated into the 
body system. The simplest choice for B is body-to-equatorial if only spherical or zonals are considered, and body-to-earth- 
fixed if tesserals are considered. 
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The operation, [p x ] . may be considered a matrix, where 

0 -p. +p t 


fr x ] s 


P* r y] 
+ Pz 0 -P* 

+ p x 0 


L~Py 


(8-27) 


then 


px G = 


0 

+ P* 
-Pu 


-p* 

0 

+ p. 


+ PL 

-p, 

0 


y 

9n 9i2 9i3 



921 922 923 

= 


931 932 ^33 



03lPy _ 02lP* 032 Py - 022 P Z 033Py _ 023P z | 
01lP z ~ 03lP z 012P Z " 032P* 9xzP x -933P 
l9 2 lP X -9nPy 9 22 P X ~ 9\ 2 P y 9 2 3P X ~ 013Py| 


(8-28) 


and finally 


031 PyP* 9 2 lP z P x + 93 2 Py 9 22 P z P y + 933PyPz 9 2 3P z 
0i 1 PzP x -9 3i Pi + 0i 2 P z P y - S32p*py + Ml - 0330*0. 


pxGp = 

Since the G matrix is always symmetric! 2 ^, we may write 


\9 2 iP 2 x -9uP y P x + 9 22 P jfiy 9i 2 Py + 023P*P Z “ SiaPypy 


(8-29) 


x = Jp x Gpdm = 


023 (^zz lyy) ^ fll3^«y 9i 2 I xz + Iyz(9s3 9 22 ) 
9l3 (*«” ^ ZZ ) ~9 2 3^xy + Pl2^y* + (011 “ 033^ 

9 12 (lyy ~ (**) + 9 2 3^ X z ~ 9l3?yz + (xy (022 ” 01 1 ) 


(8-30) 


8.3 Formulation Validation 

As a check on the general form given by eq(8-30). consider the gravity vector (assuming 
spherical planet) given by 


and consequently. 


Therefore. 


9 = -P -3 

r 


p = _ii /+ 3^ 

dr r* r 5 


G = B^B = - ^ + 3^B T n T B 
3r r* r 5 


but we defined f> to be equal to B T f. Therefore, 


G.-^,3> T 
r r 


or. 


(8-31) 


(8-32) 


(8-33) 


(8-34) 
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(8-35) 


bj - - bj b 2 bj b 3 
bj b 2 b^ - g ^ 3^2 

b,b 3 b 2 b 3 bf-i 

Substituting the G elements from eq(8-35) into eq(8-30) gives 

3)i. ^ ^3 ^ zz ~ + ^ — ^1 ^2^xz 

bjbg ( I** — I zz ) —b 2 b 3 I x y +bjb 2 /y Z 

bi b 2 (^yy — + ^2 ^3^xz ^ 3 Iy Z 

which is identical to eq(8-16). This lends confidence that the general expression is cor- 
rect and that the spherical case is contained in the more general expansion as a special 
case where all higher gravitational harmonics are zero. 

The formulae given here were coded in ASDS [9] and checked against a FORTRAN pro- 
gram containing Roithmayr's method for n=2, m=0 (J2 only). The data used and the 
results obtained are given in Appendix B. 

Further data was run for a 4x4 gravity model. Note that the gravity gradient does change 
as the model size increases. This result is expected to be especially important during 
control when longitude effects will vaiy much more quickly than those due to latitude. 


+ Iy,(bl-b*) 

+ '«<b?- b!) 
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gifl Geomagnetic Field 

The magnetic field for the earth is also defined in terms of Legendre functions. It was 
felt that it would be useful to include a derivation for the magnetic field in this report, 
since so much of the gravity algorithm can be used. 

The potential function, V, for the magnetic field is given [101, [1 11 as 

* n w+ j 

V=a^£(°) P^(cos0) (fi£'cosmX+ h^sinmX) (9-1) 

n - 1 o 

where a is the mean radius of the earth, r is position magnitude, 0 is the geocentric 
colatitude, g™ and h” are the spherical harmonic coefficients, and X is the longitude. 
The PJJ* in eq(9-l) are not derived Legendre functions, rather they are Schmidt normal- 
ized associated Legendre functions of degree n and order m. defined by 

^ (CM 0 , '[ ? ^ns!r' 5 »">] 2p ^ <c<>s6) (9 - 2) 

The cosine of the colatitude is the same as the sine of the latitude, which was denoted by 
e in eq(3-3). Also, define 


^n, m 


s 


n, m 


p2(n-m)! 

L (n+m)! 

r2 (n- m)\ 
L (n+m)! 



It is now possible to write eq(9-l) as 


(9-3) 


• n 2 n 

V= X Xt 1 ?) i, nn,(e)(C n>m cosmX + S n>m sinmX) (9-4) 

n * 1 m* 0 

This is not quite the form of the gravitational potential since the sum on n started at 2 in 
eq(3-5). Note also the appearance of a 2 in eq(9-4) in lieu of p in eq(3-5). 

Using the definitions in eq(3-6) and eq(3-9), and separating out the n = 1 term, eq(9-4) 
can be written 


up™ 


v= ^( i ? c i. 0 + (Ci.iy + s lfI y)p}) + J St ( 7 ) ^ 


(9-5) 


n * 2m* 0 


Where now, the are derived Legendre functions. The double sum part of eq(9-5) is 
now identical in form to the double sum part of eq(3-12). the only difference being that 
a 2 replaces p. Using the definitions of and P{ in eq(3-13), the leading term of eq(9-5) 
can be written 


V, = ^(P?C lt0 +(C lil - i i + S 1 . 1 ^)P}) 


£ 

r 3 


(Cj, ,X, + S,. 1 JSj + Cj flXg) 


The magnetic potential function can now be expressed as 


(9-6) 


January 1993 


NASA CR- 188243 


25 




The magnetic field vector, B. Is computed as 



(9-8) 


The partlals of the double sum part of V will look exactly as they did for gravity (assum- 
ing a 2 replaces p). That being the case, the definitions of T n , H n , J n . and K n are the 
same as given in eq(4-l 1). This allows ^ to be written 


3V = _ oV y 

BX BX r 2 l A 



(9-9) 


Computing the partial of Vj results in 


_ a* 
BX r 5 


I 

*i.i 

L.C 1 . 0 J 


-3— (Cj jXj + S 1(1 X2+ C 1 


(9-10) 


The negative of this result can be found in [121 as the field resulting from a magnetic 
dipole. Now, defining 


Jb 7 c m + X< 7>" j « 

n«2 
n * 2 

r*i (3Cj jXj + 3S 1 jX, + 2Cj 0 j%) + ^ (— ) T n 
’ <1 

”=- r Cuo+t(-/ H n 


n-2 


dV 

The equation for — can be written 
o2£ 


n * 2 


av 

a# 



r 2 


He* + ^ 


J 

K 

H 


Finally, then, the equation for the magnetic field vector, B. becomes 



r 2 


(r +zH)X- 

ji- 

ff 


_H]J 



J 

K 

Ffj 


(9-11) 


(9-12) 


(9-13) 
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Note that eq(9-13) is the same as eq(4-13) with -\i replaced by a 2 . Assuming that equal 
degree and order were desired, a great deal of time could be saved by computing the 
gravitational acceleration and the magnetic field vector together, since they both would 
use the same state, the same derived Legendre functions, and all of the sums would be 
of exactly the same form. An examination of the code given in Appendix A for the 
fast_gravity_model and the fast_magnetic_model will show that they are almost entirely 
the same. The magnetic field vector resulting from IGRF 1985 data at a given position 
vector is given in Appendix B. 
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10,0 Conclusions 


Derivations of the first and second partials of the gravitational potential have been given 
In both normalized and unnormalized form. Two different recursion relationships for the 
derived Legendre functions were considered. Code for both a “fast" unnormalized gravity 
model and a normalized gravity model using the best recursion relationship were given 
In Appendix A. Speed comparisons made in [13] indicate that the model in [2) ran at 
essentially the same speed as other gravity models in general use at JSC. The plots In 
Appendix B show that the “fast" algorithm Is always faster than the model In [2], In some 
cases by as much as 20%, and consequently should be that much faster than other 
models In general use at JSC. The normalized models are faster than the model In [2] out 
through degree and order 27. Beyond that, the extra multiplications inherent In the nor- 
malized approach begin to outweigh the savings gained by precomputing all possible 
terms. Gravitational acceleration computed by the normalized and unnormallzed mod- 
els agreed through the 15th (out of 16) significant digit, out through degree and order 
50, for a variety of Initial states. For larger models, one would probably be safer using 
code based on the recursion relationship from [3]. 

A gravity gradient torque derivation for a general gravity model was given as well as code 
and data showing that the torque agrees with another model restricted to J2 only. 

A derivation of the magnetic field vector was given. The derivation was given In a form 
that was as close as possible to the form of the gravity derivation. An examination of the 
resulting code shows that the two algorithms are almost totally alike. 
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Appendix A 


Ada Code 
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A.1 Spec of Fa*t_Gravity_Model 

with Real_Types; 
use Real_Types; 
with Vector_Matrlx_3; 
use Vector_Matrlx_3; 

package fast_Gravlty_Model Is 

Max_Gravlty_Model_N ame_Length : constant Positive := 80; 
max_degree_and_order : constant Positive := 50; 

type Data_Coefficient_Array Is 
array (Natural range <>, Natural range <>) of Real; 

type gravity_array Is array(0. .max_degree_and_order+2) of real; 
type ga_ptr is access gravity_array; 

type gravlty_array_2 is array(0..max_degree_and_order) of ga_ptr; 
type Gravlty_Model_Data Is private; 


function Create_Gravlty_Model (Name_In : String; 

C, S : Data_CoefRcient_Array; 

Mu, Radius : Real) return Gravlty_Model_Data; 


procedure Gotpot (Gmd : In Gravtty_Model_Data; 

X : In Vector_3; 

R : In Real; 

Want_Central_Force : In Boolean; 

Nax, Max : In Natural; 

G : out Vector_3); -no potential 


procedure Gotpot (Gmd : In Gravlty_Model_Data; 

X : In Vector_3; 

R : In Real; 

Want_Central_Force : In Boolean; 

Nax, Max : In Natural; 

Pot : out Real; 

G : out Vector_3; 

Dgdx : out Matrix_3x3);~ pot and dgdx 

private 

type Gravity_Model_Data Is — defaulted to point mass gem_9 
record 

Name : String (1 .. Max_Gravity_Model_Name_Length); 

Name_Length : Integer, 

C : gravity_array_2; 

S : gravlty_array_2; 

Mu : Real := 398_600.47E9; — planet gravitational constant(m**3/s**2) 

Radius : Real := 6_378_139.0; -- planet equatorial radius (m) 
Model_Max_Slze :Natural; -- max size current gravity model data 

end record; 

end fast_Gravity_Model; 
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A.2 Body of Fast_Gravlty_Model 

with Extended J*ange_Comblnatoric_Functions; 
use Extended_Range_Combinatoric_Functlons; 
with Exponentlal_Logarlthm_Functlons: 
use Exponential_Logarlthm_Functions; 

package body fast_Gravlty_Model Is 

Default_Gmd : Gravity_Model_Data; 

Have_Set_Default_Gravlty : Boolean := False; 
Gravlty_Model_Name_Too_Long : exception; 
bad_gravlty_data : exception; 
twonm 1 .twonm 1 on.nm 1 on : gravity_array; 

P : gravity_array_2 := (others => new gravity. array); 


procedure Gotpot (Gmd : In Gravity_Model_Data; 
X : in Vector_3; 

R : In Real; 

Want_Central_Force : In Boolean; 

Nax, Max : In Natural; 

G : out Vector_3) is 


Ri, Xovr. Yovr, Zovr, Ep : Real; 

Muor, Muor2, Reor, Reom,Sum_Inlt : Real; 
ctll, stll : gravlty.array; 

Sumh, Sumgam, SumJ. Sumk, Sumh.N, Lambda : Real; 
pnm.cnm, snm.ctmrn 1 ,stmm 1 : real; 

Sumgam.N, Sumj.N, Sumk.N, Mxpnm, Npmpl : Real; 
Bnmtil ,n_const : Real; 

Mml, Mm2, Mpl, Nml, Llm ,nm2: Integer; 
pn,pnml,pnm2 : ga_ptr; 
cn.sn : ga_ptr; 

begin 

Ri := 1.0 / R; 

Xovr := X (1) * RI; 

Yovr := X (2) • RI; 

Zovr := X (3) * RI; 

Ep := Zovr; 

Reor ;= gmd. Radius * Ri; 

Reom := Reor; 

Muor := gmd.Mu * RI; 

Muor2 := Muor * Ri; 

Case Want_Central_Force is 
When true => Sum.Init := 1.0; 

When false => Sum.Inlt := 0.0; 
end case; 


ctll (0) ;= 1.0; ctll (1) := Xovr; 
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stll (0) := 0.0; stll (1) := Yovr; 

Sumh := 0.0: 

Sumj := 0.0; 

Sumk := 0.0; 

Sumgam := Sum_Inlt; 

p(l)(0) := ep; 
for N in 2 .. Nax loop 
n_const := twonml(n); 

Reom := Reom * Reor; 
pn := p(n); 
cn := gmd.c(n); 
sn ;= gmd.s(n); 
rani := n - 1; 
nm2 := n - 2; 
pnml := p(nml); 
pnm2 := p(nm2); 

Pn(nml) := ep*Pn(n); 

Pn(0) := Twonmlon(n)*Ep*Pnml(0) - Nmlon(n)*Pnm2(0); 

Pn(l) := Pnm2(l) + n_const * Pnml(O); 

Sumh.N := Pn (1) * Cn(0); 

Sumgam_N := Pn (0) * Cn(0) * real(n + 1); 

if Max > 0 then 
for m In 2„nm2 loop 

Pn(m) := Pnm2(m) + n_const • Pnml(m-1): 
end loop; —Have all derived Legendre functions 
SumJ_N := 0.0; 

Sumk_N := 0.0; 

ctil (N) := ctil (1) * ctil (Nml) - stil (1) * stil (Nml); 
stil (N) := stil (1) * ctil (Nml) + ctil (1) * stil (Nml); 

if N < Max then 
Lim := N; 
else 

Lim := Max; 
end if; 

for M in 1 .. Lim loop 
Mml := M - 1; 

Mpl := M + 1; 

Npmpl := Real (N + Mpl); 

pnm := pn(m); 

cnm := cn(m); 

snm := sn(m); 

ctmml := ctil(mml); 

stmml := stil(mml); 

Mxpnm := Real (M) • Pnm; 

Bnmtil := Cnm * ctil (M) + Snm * stll (M); 

Sumh_N := Sumh_N + Pn(mpl) * Bnmtil; 

Sumgam_N := Sumgam_N + Npmpl * Pnm * Bnmtil; 

Sumj_N := SumJ_N + Mxpnm * (Cnm*ctmml + Snm*stmml); 
Sumk_N := Sumk_N - Mxpnm * (Cnm*stmml - Snm*ctmml): 
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ri r; 


end loop: 

Sumj := SumJ + Reom * SumJ_N; 

Sumk := Sumk + Reom * Sumk_N; 
end if; 

— SUMS BELOW HERE HAVE VALUES WHEN M := 0 

Sumh := Sumh + Reom * Sumh_N; 

Sumgam := Sumgam + Reom * Sumgam_N; 
end loop: 

Lambda := Sumgam + Ep • Sumh: 

G (1) := -Muor2 * (Lambda * Xovr - Sumj): 

G (2) := -Muor2 * (Lambda * Yovr - Sumk); 

G (3) := -Muor2 * (Lambda * Zovr - Sumh); 

end Gotpot: 


procedure Gotpot (Gmd : In Gravlty_Model_Data; 

X : In Vector_3; 

R : In Real; 

Want_Central_Force : in Boolean; 

Nax, Max : In Natural; 

Pot : out Real; 

G : out Vector_3: 

Dgdx : out Matrlx_3x3) is 

Mu : Real renames Gmd.Mu: 

Radius : Real renames Gmd. Radius; 

Ri, Xovr, Yovr. Zovr. Ep.Sum_Init.n_const : Real; 

Muor, Muor2. Muci3. Reor. Reom. Sumv. Gg, Ff. Dl. D2 : Real: 
ctil. stll : gravlty_array; 

Sumh, Sumgam, Sumj, Sumk, Sumh_N. Npl, Lambda : Real; 

Suml, Summ, Sumn, Sumo. Sump. Sumq, Sumr, Sums. Sumt : Real; 
Suml_N, Summ.N, Sumn_N, Sumo_N. Sump_N, Sumq_N ; Real; 
Sumr_N, Sums_N, Sumt_N, Temp : Real; 

Sumgam_N, Sumj_N, Sumk_N, Mxpnm, Npmpl : Real; 

Sumv_N. Amntll, Bnmtll, Pnmbnm, Anmtml, Bnmtml : Real; 

Mml, Mm2. Mpl, Mp2. Nml, Nm2 , Llm : Integer; 
pnm . pnmpl, cnm.snm, ctmml, stmml.cnO : real; 
pn.pnml,pnm2,cn,sn : ga_ptr; 

begin 

Rl := 1.0 / R; 

Xovr :=X (1) • Rl; 

Yovr := X (2) * Rl; 

Zovr := X (3) * Rl; 

Ep := Zovr; 

Reor := Radius * Ri; 

Reom := Reor; 

Muor := Mu * Rl; 
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Muor2 := Muor • Ri; 

Case Want_Centxal_Force Is 
When true => Sumjnlt := 1.0; 
When false => Sumjnlt := 0.0; 


end case; 

ctil (0) 

:= 1.0; ctil (1) := Xovr 

stil (0) 

:= 0.0; stil(l) := Yovr: 

Sumv 

:= Sum Inlt; 

Sumh 

:= 0.0; 

Sum] 

:= 0.0; 

Sumk 

:= 0.0; 

Sumgam ;= Sumjnlt; 

Summ 

:= 0.0; 

Sumn 

:= 0.0; 

Sumo 

:= 0.0; 

Sump 

:= 0.0; 

Sumq 

:= 0.0; 

Sumr 

:= 0.0; 

Sums 

:= 0.0; 

Sumt 

:= 0.0; 

Suml 

:= 2.0 * Sumjnlt; 

P(1)(0) 

:= ep; 

for N in 2 .. Nax loop 


n_const ;= Twonml(n); 

Reom := Reom • Reor; 
nml := n - 1; 
nm2 := n - 2; 
pn := p(nj; 
pnml := p(nml); 
pnm2 := p(nm2); 

Pn(nml) := ep*Pn(n); 

Pn(0) := Twonmlon(n)*Ep*Pnml(0) - Nmlon(n)*Pnm2(0); 
Pn(l) := Pnm2(l) + n_const * Pnml(0); 
cn := gmd.c(n); 
sn := gmd.s(n); 
npl := real(n+l); 

CnO := Cn(0); 

Sumv_N ;= Pn (0) * CnO; 

Sumh_N := Pn (1) • CnO; 

Summ_N := Pn (2) * CnO; 

Sumgam_N := Sumv_N * Npl; 

Sump_N := Sumh_N • Npl; 

Suml_N := Sumgam_N • (Npl + 1.0); 

If Max > 0 then 
for m in 2..nm2 loop 

Pn(m) := Pnm2(m) + n_const * Pnml(m-l); 
end loop; 

Sumj_N := 0.0; 

Sumk_N := 0.0; 
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Sumn_N := 0.0; 

Sumo_N ;= 0.0; 

Sumq_N ;= 0.0; 

Sumr_N := 0.0; 

Sums_N := 0.0; 

Sumt_N := 0.0; 

nml := n - 1; 

ctll (N) := ctll (1) * ctil (Nml) - stll (1) * stil (Nml); 

stll (N) := stil (1) * ctll (Nml) + ctll (1) * stll (Nml); 

If N < Max then 
Llm := N; 
else 

Llm := Max; 
end If; 

for M In 1 .. Llm loop 
Mml := M - 1 ; 

Mpl := M + 1; 

Mp2 := M + 2; 

Npmpl := Real (N + Mpl); 

pnm := pn(m); 
pnmpl := pn(mpl); 
cnm := cn(m); 
snm := sn(m); 
ctmml := ctil(mml); 
stmml := stll(mrol); 

Mxpnm := Real (M) • pnm; 

Bnmtll := cnm * ctll (M) + snm * stll (M); 

Pnmbnm := Pnm * Bnmtll; 

Sumv_N := Sumv_N + Pnmbnm: 

Bnmtml := CNm * ctMml + SNm * stMml; 

Anmtml := CNm * stMml - SNm * ctMml; 

Sumh_N := Sumh_N + Pn (Mpl) * Bnmtll; 

Sumgam_N := Sumgam_N + Npmpl * Pnmbnm; 

SumJ_N := SumJ_N + Mxpnm * Bnmtml; 

Sumk_N := Sumk_N - Mxpnm * Anmtml; 

Suml_N ;= Suml_N + Npmpl * (Real (Mpl) + Npl) ‘pnmBnm; 
Summ_N := Summ_N + Pn(Mp2) * Bnmtll; 

Sump_N ;= Sump_N + Npmpl * PnMpl * Bnmtll; 

Sumq_N := Sumq_N + Real (M) * PnMpl * Bnmtml; 

Sumr_N := Sumr_N - Real (M) * PnMpl * Anmtml; 

Sums_N ;= Sums_N + Npmpl * Mxpnm * Bnmtml; 

Sumt_N := Sumt_N - Npmpl * Mxpnm * Anmtml; 

If (M >= 2) then 
Mm2 := M - 2; 

Sumn_N := Sumn_N + Real (Mml) * Mxpnm * 

(CNm * ctll (Mm2) + SNm*stll(Mm2)); 

Sumo_N := Sumo_N + Real (Mml) * Mxpnm • 

(CNm • stll (Mm2) - SNm*ctil{Mm2)): 

end if; 
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end loop: 

Sumj := SumJ + Reom * SumJ_N; 

Sumk := Sumk + Reom * Sumk_N ; 

Sumn := Sumn + Reom * Sumn_N : 

Sumo := Sumo + Reom * Sumo_N; 

Sumq := Sumq + Reom * Sumq_N; 

Sumr := Sumr + Reom * Sumr_N; 

Sums := Sums + Reom * Sums_N; 

Sumt := Sumt + Reom * Sumt_N; 
end If: 

— SUMS BELOW HERE HAVE VALUES WHEN M := 0 
Sumv := Sumv + Reom * Sumv_N; 

Sumh := Sumh + Reom * Sumh_N; 

Sumgam := Sumgam + Reom * Sumgam_N; 

Sum] := Sum! + Reom * Suml_N; 

Summ := Summ + Reom * Summ_N; 

Sump := Sump + Reom * Sump_N; 
end loop; 

Pot := Muor • Sumv; 

Lambda := Sumgam + Ep * Sumh; 

G (1) := -Muor2 * (Lambda * Xovr - Sumj); 

G (2) := -Muor2 • (Lambda * Yovr - Sumk); 

G (3) := -Muor2 • (Lambda * Zovr • Sumh); 

-- Need to construct second partial matrix_3x3 
Gg := -(Summ * Ep + Sump + Sumh); 

Ff := Sum] + Lambda + Ep * (Sump + Sumh - Gg): 

D 1 := Ep * Sumq + Sums; 

D2 := Ep • Sumr + Sumt; 

Muor3 := Muor2 • Ri; 

Dgdx (1. 1) := Muor3 * ((FT* Xovr - 2.0 • Dl) * Xovr - Lambda + Sumn); 
Dgdx (2. 2) := Muor3 * ((Ff • Yovr - 2.0 * D2) • Yovr - Lambda - Sumn); 
Dgdx (3, 3) := Muor3 • ((Ff • Zovr + 2.0 * Gg) • Zovr - Lambda + Summ); 
Temp ;= Muor3 * ((Ff • Yovr - D2) * Xovr - Dl * Yovr - Sumo); 

Dgdx (1,2):= Temp; 

Dgdx (2, 1) := Temp; 

Temp := Muor3 * ((Ff • Xovr - Dl) * Zovr + Gg • Xovr + Sumq); 

Dgdx (1,3):= Temp; 

Dgdx (3. 1) := Temp; 

Temp := Muor3 * ((Ff * Yovr - D2) * Zovr + Gg * Yovr + Sumr); 

Dgdx (2. 3) := Temp; 

Dgdx (3, 2) ;= Temp; 

end Gotpot; 


function Create_Gravlty_Model (Name_In : String: 

C. S : Data_Coefficient_Array; 

Mu, Radius : Real) return Gravity _Model_Data is 
Gmd : Gravity_Model_Data; 

Coef : Real; 

n_max : Integer := CLast (1); 
begin 
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If n_max < 2 then raise bad_gravity_data;end if; 

gmd.c := (others => new gravity_array); 
gmd.s := (others => new gravity_array); 

— Unnormalize gravity model coefficients 
for N in CRange loop 
for M in 0 .. N loop 
if M = 0 then 

Gmd.C (N)(0) := SqrtfReal (2 * N + 1)) * C (N. 0) * 1.0E-6; 
Gmd.S (N)(0) := 0.0; 
else 

Coef := 

Sqrt (Real (2 * (2 * N + 1)) * Factorial_Ratio (N - M. N + M)) * 
1.0E-6; 

Gmd.C (N)( M) := Coef • C (N, M); 

Gmd.S (N)( M) := Coef * S (N. M); 
end if; 
end loop; 
end loop; 

Gmd.Mu ;= Mu; 

Gmd.Radius := Radius; 

Gmd.Name_Length ;= Name_In’Length; 

if Gmd.Name_Length > Max_Gravity_Model_Name_Length then 
raise Gravity_Model_Name_Too_Long; 
end if; 

Gmd.Name := (others => Ascil.Nul): 

Gmd.Name (1 .. Gmd.Name.Length) := Name_In; 
Gmd.Model_Max_Size := njnax; 

if Have_Set_Default_Gravity then 
return Gmd; 
else 

Default_Gmd := Gmd; 
return Gmd; 
end if; 

end Create_Gravity_Model; 


begin 


p(0)(0) ;= 1.0; p(0)(l) := 0.0; p(0)(2) := 0.0; 
p(l)(l) := 1.0; p(l)(2) := 0.0; p(l)(3) := 0.0; 
for n in 2..Max_Degree_And_Order loop 
p(n)(n) := p(n-l)(n-l)*real(2*n-l); 
p(n)(n+l) := 0.0; 
p(n)(n+2) := 0.0: 
twonml(n) := real(2*n - 1); 
twonmlon(n) := twonml(n)/real(n); 
nmlon(n) := real(n - l)/real(n); 
end loop: 

end fast_Gravity_Model; 
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A.3 Body of Norma llzed_Gravity_Model 

with mathematical_constants; 
with exponential_logarlthm_functlons; 
use exponential_logarlthm_functions; 

package body Normalized_Gravity_Model is 

sqrt3:constant real := mathematlcal_constants.square_root_of_three; 
Gravlty_Model_Name_Too_Long : exception; 
bad_gravity_data : exception; 

P : gravity_array_2 := (others => new gravlty_array); 
pi :ga_ptr:= p(l); 

xl : gravtty_array_2 := (others => new gravity_array); 
eta : gravlty_array_2 := (others => new gravlty_array); 
zeta : gravlty_array_2 := (others => new gravity_array); 
ups lion: gravity_array_2 := (others => new gravlty_array); 
alpha : gravlty_array ; 
beta : gravlty_array ; 
nrdlag : gravlty_array ; 
num.den: Integer; 


function Create_Gravlty_Model (Name_In : String; 

C. S : D ata_Coefficient_Array; 

Mu. Radius : Real) return Gravity_Model_Data Is 
Gmd : Gravlty_Model_Data; 

Coef : Real; 

n_max : Integer := CXast (1); 
begin 

If n_max < 2 then raise bad_gravlty_data;end If; 

gmd.c := (others => new gravity_array); 
gmd.s := (others => new gravlty_array); 

— scale gravity model coefficients 
for N In C ’Range loop 
for M In 0 .. N loop 

Gmd.C (N){ M) := 1.0e-6 • C (N. M);~Just scale coefficients 
Gmd.S (N)( M) := I.0e-6 * S (N. M);-- Just scale coefficients 
end loop; 
end loop; 

Gmd. Mu := Mu; 

Gmd.Radlus := Radius; 

Gmd.NameJLength := NameJn’Length; 

If Gmd.Name_Length > Max_Gravity_Model_Name_Length then 
raise Gravlty_Model_Name_Too_Long; 
end if; 

Gmd.Name := (others => Ascii.Nul); 

Gmd. Name (1 .. Gmd.Name_Length) := Name_In; 
Gmd.Model_Max_Slze := n_max; 

return Gmd; 
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end Create Jjravity_Model; 


procedure Gotpot (Gmd : In Gravlty_Model_Data; 
X : in Vector_3; 

R : In Real; 

Want_Central_Force : In Boolean; 

Nax, Max : In Natural; 

G ; out Vector_3) is 


Ctil, Stil ; gravity_array; 

Ri, Xovr, Yovr, Zovr, Ep, Sumjnit : Real; 

Muor, Muor2, Reor, Reom : Real; 

Sumh, Sumgam, SumJ, Sumk, Npl, Lambda ; Real; 

Sumh_N,Sumgam_N, Sumj_N, Sumk_N, Mxpnm, Npmpl : Real; 

Bnmtil.pnm.snm.cnm.ctmml.stmml : Real; 

Mini, Mpl, Nml, nm2,Lim: Integer; 

pn,pnml,pnm2 : ga_ptr; 

cn,sn,zn^dn,etn : ga_ptr; 

begin 

Ri := 1.0 / R; 

Xovr := X (1) * Ri; 

Yovr ;= X (2) * Ri; 

Zovr := X (3) * Ri; 

Ep := Zovr; 

Reor ;= gmd. Radius * Ri; 

Reom := Reor; 

Muor := gmd.Mu * Ri; 

Muor2 := Muor * Ri; 

Case Want_Central_Force is 
When true => Sumjnit := 1.0; 

When false => Sumjnit := 0.0; 

end case; 

Ctil (0) ;= 1.0; Ctil (1) := Xovr. 

Stil (0) := 0.0; Stil (1) := Yovr; 

Sumh := 0.0; 

SumJ := 0.0; 

Sumk := 0.0; 

Sumgam := Sumjnit; 

pl(0) := sqrt3*ep; 

for N in 2 .. Nax loop 
Reom ;= Reom • Reor; 
pn := p(n); 
cn := gmd.cfn); 
sn ;= gmd.s(n); 
zn := zeta(n); 
xin ;= xi(n); 
etn ;= eta(n); 
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nml := n-1; 
nm2 := n-2; 
pnml := p(nml); 
pnm2 := p(nm2); 

Pn(0) := alpha(n)*Ep*Pnml(0) - beta(n)*Pnm2(0); 
Pn(nml) := ep*nrdlag(n); 

Pn(l) :■ xln(l)*ep*Pnml(l) - etn(l)* Pnm2(l): 
Sumh_N := zn(0)*Pn(l) * Cn(0); 

Sumgam_N := Pn(0) * Cn(0) * real[n + 1); 

if Max > 0 then 
for m in 2..nm2 loop 

Pn(m) := xtn(m)*ep*Pnml(m) - etn(m)* Pnm2(m); 
end loop: --got all the Legendre functions now 

Sumj_N := 0.0; 

SumK_N := 0.0; 

Ctll (N) := Ctll (1) * Ctil (Nml) - Stll (1) * Stil (Nml); 
Stll (N) := Stll (1) * Ctll (Nml) + Ctll (1) * Stll (Nml); 


If N < Max then 
Lim := n; 
else 

lim := Max; 
end If; 

for M In 1 .. Llm loop 
Mml := M - 1; 

Mpl := M + 1; 

Npmpl := Real (N + Mpl); 

pnm := pn(m); 

cnm := cn(m); 

snm := sn(m); 

ctmml := ctil(mml); 

stmml := stllfmml); 

Mxpnm := Real (m) * Pnm; 

Bnmtll := Cnm * Ctll (M) + Snm * Stll (M); 

Sumh_N := Sumh_N + Pn(mpl) * Bnmtll*zn(m); 

Sumgam_N := Sumgam_N + Npmp 1 • Pnm * Bnmtll; 

Sumj_N := Sumj_N + Mxpnm*(Cnm*ctmm 1 + Snm*stmml); 
Sumk_N := Sumk_N - Mxpnm*(Cnm*stmml - Snm*ctmml); 
end loop: 

SumJ := SumJ + Reom * Sumj_N; 

Stunk := Sumk + Reom * Sumk_N; 
end If; 

— - SUMS BELOW HERE HAVE VALUES WHEN M := 0 

Sumh := Sumh + Reom * Sumh_N; 

Sumgam ;= Sumgam + Reom * Sumgam_N; 

end loop; 
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in: 


Lambda := Suing am + Ep • Sumh; 

G (1) := -Muor2 * (Lambda * Xovr - Sumj); 
G (2) := -Muor2 * (Lambda * Yovr - Sumk); 
G (3) := -Muoi2 * (Lambda * Zovr - Sumh); 

end Gotpot; 


procedure Gotpot (Gmd : In Gravlty_Model_Data; 

X : In Vector_3; 

R : In Real; 

Want_C e ntral_Forc e : in Boolean; 

Nax, Max : in Natural; 

Pot ; out Real; 

G ; out Vector_3; 

Dgdx : out Matrix_3x3) Is —pot and dgdx 
Ctll, Stll : gravlty_array; 

Rl. Xovr, Yovr, Zovr, Ep,Sum_Inlt ; Real; 

Muor, Muor2, Muor3, Reor, Reom, Sumv, Gg, Ff, Dl, D2 : Real; 
Sumh, Surngam. Sumj, Sumk, Sumh_N, Npl, Lambda : Real: 

Suml. Summ, Sumn, Sumo, Sump, Sumq, Sumr, Sums, Sumt : Real; 
Suml_N. Summ_N, Sumn_N, Sumo_N, Sump_N, Sumq_N : Real; 
Sumr_N. Sums_N, Sumt_N, Temp : Real; 

Sumgam_N, SumJ_N, Sumk_N, Mxpnm, Npmpl : Real; 

Sumv N, Amntll, Bnmtil, Pnmbnm, Anmtml, Bnmtml : Real; 

Mini, Mm2, Mpl, Mp2, Nml, Nm2 , Llm : Integer; 
pnm.pnmpl.cnm.snm, ctmml, stmml,z_pnmpl,cnO : real; 
pn.pnml.pnm2 : gajptr; 
cn,sn.zn,upsn,xln,etn : ga_ptr; 


begin 

Rl := 1.0 / R; 

Xovr := X (1) * Rl; 

Yovr := X (2) * Rl; 

Zovr := X (3) * Rl; 

Ep ;= Zovr, 

Reor := gmd. Radius * Rl; 

Reom := Reor: 

Muor := gmd.Mu * Rl; 

Muor2 := Muor * Rl; 

Muor3 := Muor2 * Rl; 

Case Want_Central_Force Is 
When true => Sum_Init ;= 1.0; 
When false => Sum_Init := 0.0; 
end case; 

ctll (0) := 1.0; ctll (1) := Xovr; 
stll (0) := 0.0; stll (1) := Yovr; 
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Sumv 

:= Sum_Inlt; 

Sumh 

:= 0.0; 

SumJ : 

= 0.0; 

Sumk 

:= 0.0; 

Sumgam 

:= Sum Init; 

Summ 

:= 0.0; 

Suran 

:= 0.0; 

Sumo 

:= 0.0; 

Sump 

:= 0.0; 

Sumq 

:= 0.0; 

Sumr 

:« 0.0; 

Sums 

:= 0.0; 

Sumt 

= 0.0; 

Sum! ; 

= 2.0 * Sum_Inlt; 

P(1)(0) := 

sqrt3*ep: 

for N In 2 

.. Nax loop 


Reom := Reom • Reor; 

pn := p(n): 

cn := gmd.c(n); 

sn := gmd.s(n); 

zn := zeta(n): 

xin := xl(n); 

etn := eta(n): 

nml := n - 1; 

mn2 := n - 2; 

pnml := p(nml): 

pnm2 := p(nm2); 

Pn(0) := alpha(n) *Ep*Pnm 1 (0) - beta(n)*Pnm2(0); 
Pn(nml) := ep*nrdJag(n); 

Pn(l) := xin(l)*ep*Pnml(l) - etn(l)* Pnm2(l); 
upsn := upsilon(n); 
npl := real(n+l); 

CnO := Cn(0); 

Sumv_N := Pn (0) * CnO; 

Sumh_N := Pn (1) * CnO *zn(0); 

Summ_N := Pn (2) * Cn0*upsn(0); 

Sumgam_N := Sumv_N * Npl; 

Sump_N := Sumh_N * Npl; 

Suml_N := Sumgam_N • (Npl + 1.0): 

If Max > 0 then 
for m In 2..nm2 loop 

Pn(m) := xln(m) *ep*Pnm 1 (m) - etn(m)* Pnm2(m); 
end loop; —got all the Legendre functions now 

Sumj_N := 0.0; 

Surnk^N := 0.0; 

Sumn_N ;= 0.0; 

Sumo_N := 0.0; 

Sumq_N := 0.0; 

Sumr_N := 0.0; 

Sums_N ;= 0.0; 
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Sumt_N := 0.0; 


ctil (N) := ctil (1) * ctll (Nml) - stll (1) * stll (Nml); 
stll (N) := stll (1) * ctil (Nml) + ctil (1) * stll (Nml); 

if N < Max then 
Lim ;= N; 
else 

Lim := Max; 
end if; 

for M in 1 .. Lim loop 
Mml := M - 1; 

Mpl := M + 1; 

Mp2 := M + 2; 

Npmpl := Real (N + Mpl); 

pnm := pn(m); 
pnmpl := pn(mpl); 
cnm ;= cn(m); 
snm := sn(m); 
ctmml := ctll(mml); 
stmml := stil(mml); 

Mxpnm := Real (m) * Pnm; 

Bnmtil := Cnm * Ctil (M) + Snm * Stll (M); 

Bnmtml := Cnm • ctMml + Snm * stMml ; 

Anmtml := Cnm * stMml - Snm * ctMml; 

Pnmbnm := Pnm * Bnmtil; 

Sumv_N := Sumv_N + Pnmbnm; 

If m < n then 

z_pnmpl :=zn(m)*Pn(mpl); 

Sumh_N := Sumh_N + z_pnmpl * Bnmtil; 

Sump_N := Sump_N + Npmpl • z_PnMpl • Bnmtil; 
Sumq_N := Sumq_N + Real (M) * z_PnMpl * Bnmtml; 
Sumr_N := Sumr_N - Real (M) * z_PnMpl * Anmtml; 
end if; 

Sumgam_N := Sumgam_N + Npmpl • Pnmbnm: 

Sumj_N ;= Sumj_N + Mxpnm *bnmtm 1 ; 

Sumk_N := Sumk_N - Mxpnm'anmtml; 

Suml_N := Suml_N + Npmpl * (Real (Mpl) + Npl) ‘pnmBnm; 
Summ_N := Summ_N + Pn(Mp2) * Bnmtll*upsn(m) ; 

Sums_N := Sums_N + Npmpl * Mxpnm * Bnmtml; 

Sumt_N := Sumt_N - Npmpl ♦ Mxpnm * Anmtml; 

If (M >= 2) then 
Mm2 := M - 2; 

Sumn_N := Sumn_N + Real (Mml) * Mxpnm * 

(Cnm * ctil (Mm2) + Snm*stll(Mm2)) ; 

Sumo_N := Sumo_N + Real (Mml) * Mxpnm * 

(Cnm * stll (Mm2) - Snm*ctil(Mm2)); 

end If; 
end loop; 
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Sumj := Sumj + Reom * Sumj_N; 

Sumk := Sumk + Reom * Sumk_N; 

Sum := Sum + Reom * Sumn_N; 

Sumo := Sumo + Reom * Sumo_N; 

Sumq := Sumq + Reom * Sumq_N; 

Sumr := Sumr + Reom * Sumr_N; 

Sums := Sums + Reom * Sums_N; 

Sumt := Sumt + Reom * Sumt_N; 
end If; 

— SUMS BELOW HERE HAVE VALUES WHEN M := 0 

Sumv := Sumv + Reom • Sumv_N; 

Sumh := Sumh + Reom * Sumh_N; 

Sumgam := Sumgam + Reom * Sumgam_N; 

Suml := Suml + Reom • Suml_N; 

Summ := Summ + Reom * Summ_N; 

Sump := Sump + Reom * Sump_N; 

end loop; 

Pot := Muor • Sumv; 

Lambda := Sumgam + Ep * Sumh; 

G (1) := -Muor2 * (Lambda * Xovr - Sumj); 

G (2) := -Muor2 * (Lambda * Yovr - Sumk); 

G (3) := -Muor2 * (Lambda * Zovr - Sumh); 

— Need to construct second partial matrix_3x3 
Gg := -(Summ * Ep + Sump + Sumh); 

Ff := Suml + Lambda + Ep * (Sump + Sumh - Gg); 

D1 := Ep * Sumq + Sums; 

D2 := Ep • Sumr + Sumt; 

Dgdx (1, 1) := Muor3 * ((Ff * Xovr - 2.0 * Dl) * Xovr - Lambda + Sumn); 

Dgdx (2, 2) := Muor3 • ((Ff * Yovr - 2.0 * D2) * Yovr - Lambda - Sumn); 

Dgdx (3. 3) := Muor3 * ((Ff * Zovr + 2.0 * Gg) * Zovr - Lambda + Summ); 

Temp := Muor3 * ((Ff * Yovr - D2) ‘ Xovr - Dl * Yovr - Sumo); 

Dgdx (1, 2) := Temp; 

Dgdx (2, 1) := Temp; 

Temp := Muor3 * ((Ff * Xovr - D 1) * Zovr + Gg * Xovr + Sumq); 

Dgdx (1, 3) := Temp; 

Dgdx (3. 1) := Temp; 

Temp := Muor3 * ((Ff * Yovr - D2) * Zovr + Gg * Yovr + Sumr); 

Dgdx (2, 3) := Temp; 

Dgdx (3. 2) := Temp; 

end Gotpot; 


BEGIN 

for n In 2..max_degree_and_order loop 
for m In 0 .. n- 1 loop 
num := (2*n-l)*(2*n+l); 
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den := (n+m)*(n-m); 
xl(n)(m) := sqrt{ real(num)/real(den)); 
end loop; 
end loop; 

for n In 2..max_degree_and_order loop 
for m In 0 .. n- 1 loop 
num := (2*n+l)*(n+m-l)*(n-m-l); 
den := (n+m)*(n-m)*(2*n-3); 

If num = 0 then 
eta(n)(m) := 0.0; 
else 

eta(n)(m) := sqrt( real(num)/real(den)); 
end if; 
end loop; 
end loop; 

for n In 2..max_degree_and_order loop 
for m In 0 .. n loop 
if m = 0 then 
num := n*(n+l); 

zeta(n)(0) := sqrt( real(num)/2.0); 
else 

num ;= (n-m)*(n+m+l); 
if num = 0 then 
zeta(n)(m) := 0.0; 
else 

zeta(n)(m) := sqrt(real(num)); 
end if; 
end if; 
end loop; 
end loop; 

for n in 2..max_degree_and_order loop 
for m in 0 .. n loop 
if m = 0 then 

num := n*(n-l)*(n+l)*(n+2); 
upsilon(n)(0) := sqrt( real(num)/2.0); 
else 

num := (n-m) *(n+m+ 1 )*(n-m- 1 ) *(n+m+2) ; 
if num = 0 then 
upsilon(n)(m) := 0.0; 
else 

upsilon(n)(m) := sqrt(real(num)); 
end if; 
end if; 
end loop; 
end loop; 

p(0)(0) := 1.0; p(0)(l) ;= 0.0; p(0)(2) := 0.0; 

p(l)(l) := sqrt3; p(l)(2) := 0.0; p(l)(3) := 0.0; 

for n in 2..Max_Degree_And_Order loop 

p(n)(n) ;= sqrt( real(2*n+ lT/real(2*n))*p(n- l)(n- 1); 
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nrdiag(n) := sqrt( real(2*n+l))*p(n-l)(n-l); 

num := (2*n+l)*(2*n-l); 

alpha(n) := sqrt(real(num))/real(n); 

num := (2*n+l); 

den := (2*n-3); 

beta(n) := sqrt(real(num)/real(den))*real(n-l)/real(n); 
end loop; 

end Normalized_Gravity_Model: 
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A.4 Body of General_Gravity_Gradient 

with trlgonometric.functions; use trlgonometric_functlons: 

package body general_gravity_gradient is 

function gravlty_gradlent_torque(mass_tensor, dgdx : matrlx_3x3; 
pltch.yaw,roll:rea]) return vector_3 Is 

torque : vector_3; 

g. b : matrlx_3x3 ; 

s 1 ,c 1 ,s2,c2,s3,c3,s2s3 t s2c3:real; 

begin 

si := sln(pitch); 
cl := cos(pltch); 
s2 := sinfyaw); 
c2 := cos (yaw); 
s3 := sln(roll); 
c3 := cos(roll): 

S2s3 := S2 * S3; 

S2c3 := S2 * C3: 

b := (( Cl • C2, -Cl * S2c3 + SI • S3. Cl * S2s3 + SI * C3). 

( S2 . C2*C3 ,-C2 * S3 ), 

(-S1 • C2, SI * S2c3 + Cl * S3.-S1 * S2s3 + Cl * C3)); 

g := b**tr * dgdx * b ; —Note: **tr results In transpose 

torque(l) := g(2,3)*(mass_tensor(3,3) - mass_tensor(2,2)) 

- g(l,3)* mass_tensor( 1 ,2) 

+ g{1.2)* mass_tensor( 1 .3) 

- mass_tensort2,3)*(g(3,3) - g(2.2)) ; 

torque(2) := g(l,3)*(mass_tensor(l.l) - mass_tensor(3,3)) 

+ g(2.3)* mass_tensor(l ,2) 

- g(1.2)* mass_tensor(2,3) 

- mass_t ensort 1 ,3) *(g( 1 . 1 ) - g(3.3)) : 

torque(3) := g(1.2)*(mass_tensor(2,2) - mass_tensor( 1,1)) 

- g(2,3)* mass_tensor(l,3) 

+ g(l,3)* mass_tensor(2,3) 

- massjensort 1 ,2)*(g(2,2) - g(l.l)) ; 

return torque ; 
end gravlty_gradient_torque; 
end general_gravity_gradient ; 
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A.5 Spec of Fast_Magnetlc_Model 


with Real_Types; 
use Realjiypes; 
with Vector_Matrlx_3; 
use Vector_M atrlx_3 ; 

package fast_Magnetic_Model Is 

Max_Magnetlc_Model_Name_Length : constant Positive := 80; 
max_degree_and_order : constant Positive := 20; 

type Data_Coefficlent_Array Is 
array (Natural range <>, Natural range <>) of Real; 
type magnetlc_array is array(0. .max_degree_and_order+2) of real; 
type mag_ptr Is access magnetlc_array; 
type magnetlc_array_2 is arrayfO. ,max_degree_and_order) of 

mag_ptr; 

type Magnetlc_Model_Data is private; 


function Create_Magnetic_Model (Name_In : String; 

G, H : Data_Coefficient_AiTay; 

Radius : Real) return Magnetic_Model_Data; 


procedure Maggot (Mmd : In Magnetic_Model Data; 
X : In Vector_3; 

R : In Real; 

Nax, Max ; In Natural; 

B : out Vector_3); 


private 

type Magnetic_Model_Data is — defaulted to point mass gem_9 
record 

Name : String (1 .. Max_magnetlc_Model_Name_Length); 

Name_Length : Integer; 

G : magnet lc_array_2 ; 

H : magnetlc_array_2; 

Radius : Real := 6_371_200.0; — planet mean radius (m) 
Model_Max_Slze : Natural;— max size current gravity model data 
end record; 


end fast_Magnetic_Model; 
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A.6 Body of Fast_Magnetic_Model 

with Extended_Range_Combinatoric_Functions; 
use Extended_Range_Combinatoric_Functions; 
with Exponentlal_Logarithm_Fu notions: 
use Exponential_Logarlthm_Functions: 

package body Fast_Magnetlc_Model Is 

Magnetlc_Model_Name_Too_Long : exception; 

bad_Magnetic_data ; exception; 

twonm 1 .twonm 1 on,nm 1 on : Magnetlc_array; 

P : Magnetlc_array_2 := (others => new Magnetlc_array); 


procedure Maggot (Mmd : In Magnetic_Model_Data; 

X : In Vector_3; 

R : In Real; 

Nax, Max : In Natural; 

B : out Vector_3) Is 


Rl, Xovr, Yovr, Zovr, Ep : Real; 

Muor, aeor2, Reor, Reom : Real; 
ctil, stll ; Magnetlc_array; 

Sumh, Sumgam, Sumj, Sumk, Sumh_N, Lambda : Real; 
pnm.cnm.snm.ctmm 1 .stxnm 1 : real; 

Sumgam_N, Sumj_N, Sumk_N, Mxpnm. Npmpl : Real; 
Bnmtll,n_const : Real; 

Mml, Mm2, Mpl, Mp2, Nml, Llm ,nm2: Integer 
pn,pnml,pnm2 : mag_ptr; 
cn,sn : mag_ptr; 

begin 

Rl := 1.0 / R; 

Xovr := X (1) * Rl; 

Yovr := X (2) * Rl; 

Zovr := X (3) • Rl; 

Ep := Zovr; 

Reor := Mmd. Radius • Rl; 

Reom := Reor; 
aeor2 := reor* reor; 

ctil (0) := 1.0; ctil (1) := Xovr; 
stll (0) := 0.0; stll (1) := Yovr; 

If Nax < 1 then 
Sumj := 0.0; 

Sumk := 0.0; 

Sumh := 0.0; 

Sumgam := 0.0; 
elsif Max > 0 then 
Sumj := reor*mmd.g( 1)( 1); 

Sumk := reor*mmd.h(l)(l); 
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Sumh := reor*mmd.g(l)(0); 

Sumgam := 3.0*(sumJ*xovr + sumk*yovr) + 2.0*(sumh*zovr); 
else — Max = 0 
SumJ := 0.0; 

Sumk := 0.0; 

Sumh := reor*mmd.g(l)(0); 

Sumgam := 2.0*sumh*zovr; —note: ep and zovr are the same thing 
end If; 

p(l)(0) := ep: 
for N In 2 .. Nax loop 
n_const := twonml(n); 
nml := n - 1; 
nm2 := n - 2; 
pn := p(n); 
pnml := p(nml); 
pnm2 := p(nm2); 

Pn(nml) := ep*Pn(n); 

Pn(0) := Twonmlon(n)*Ep*Pnml(0) - Nmlon(n)*Pnm2(0); 

Pn(l) := Pnm2(l) + n_const * Pnml(0); 

Reom := Reom • Reor; 
cn := mmd.g(n); 
sn := mmd.h(n); 

Sumh.N := Pn (1) * Cn(0); 

Sumgam_N := Pn (0) * Cn(0) * real(n +1); 

If Max > 0 then 
for m in 2..nm2 loop 

Pn(m) := Pnm2(m) + n_const * Pnml(m-l); 
end loop; 

SumJ_N := 0.0; 

Sumk_N := 0.0; 
nml := n - 1; 

ctll (N) := ctll (1) * ctll (Nml) - stil (1) * stll (Nml); 
stll (N) := stll (1) * ctll (Nml) + ctll (1) * stll (Nml); 

If N < Max then 
Urn := N; 
else 

Lim := Max; 
end If; 

for M In 1 .. Lim loop 
Mml := M - 1; 

Mpl := M + 1; 

Npmpl := Real (N + Mpl); 

pnm := pn(m); 

cnm := cn(m); 

snm := sn(m); 

ctmml := ctll(mml); 

stmml := stll(mml); 

Mxpnm := Real (M) * Pnm; 
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Bnmtll := Cnm * ctU (M) + Snm * stll (M); 

Sumh_N := Sumh_N + Pn(mpl) • Bnmtll; 

Sumgam_N := Sumgam_N + Npmpl * Pnm * Bnmtll; 

Sumj_N := Sumj_N + Mxpnm * (Cnm*ctmml + Snm*stmml); 
Sumk_N := Sumk_N - Mxpnm * (Cnm*stmml - Snm'ctmml); 
end loop; 

Sumj := Sumj + Reom * SumJ_N; 

Sumk := Sumk + Reom * Sumk_N; 
end If; 

— SUMS BELOW HERE HAVE VALUES WHEN M := 0 

Sumh := Sumh + Reom * Sumh.N; 

Sumgam := Sumgam + Reom * Suing am_N; 

end loop; 

Lambda := Sumgam + Ep * Sumh; 

B (1) := aeoi2 * (Lambda * Xovr - Sumj); 

B (2) := aeor2 * (Lambda * Yovr - Sumk); 

B (3) := aeor2 • (Lambda * Zovr - Sumh): 

end Maggot; 


function Create_Magnetlc_Model (Name_In : String; 

g. h ; Data_Coefficient_Array; 

Radius : Real) return Magnetlc_Model_Data is 
Gmd : Magnetic_Model_Data; 

Coef : Real; 

n_max : Integer := g'Last (1); 
begin 

If n_max < 2 then raise bad_magnetlc_data; end If; 

gmd.g := (others => new Magnetic_array) ; 
gmd.h := (others => new Magnetlc_array); 

— Unnormallze gravity model coefficients 
for N In g’Range(l) loop 
for M In 0 .. N loop 
If M = 0 then 

Gmd.G (N)(0) ;= G (N, 0)*1.0e-9 ; 

Gmd.H (N)(0) := 0.0; 
else 

Coef := Sqrt (2 .0*Factoiial Ratio(N - M.N + M))*1.0e-9; 

Gmd.G (N)( M) := Coef * G (N. M); 

Gmd.H (N)( M) := Coef * H (N. M); 
end If; 
end loop; 
end loop; 

Gmd.Radlus := Radius; 

Gmd. Name_Length := Name_In*Length; 

If Gmd.Name_Length > Max_Magnetlc_Model_Name_Length then 
raise Magnetic_Model_Name_Too_Long; 
end if; 
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Gmd.Name := (others => Ascii.Nul); 

Gmd.Name (1 .. Gmd.Name_Length) := Name_In; 
Gmd.Model_Max_Slze := n_max; 

return Gmd; 

end Create_Magnetlc_Model; 


begin —Initialize constant values 


p(0)(0) := 1.0; p(0)(l) := 0.0; p(0)(2) := 0.0; 

p(l)(l) ;= 1.0; p(l)(2) := 0.0; p(l)(3) := 0.0; 
for n In 2 . . Max_Degree_And_Order loop 
p(n)(n) := p(n- 1 )(n- 1 ) *real(2 *n- 1 ) ; 
p(n)(n+l) := 0.0; 
p(n)(n+2) := 0.0; 
twonml(n) := real(2*n - 1); 
twonmlon(n) := twonml(n)/real(n); 
nmlon(n) ;= real(n - l)/real(n); 
end loop; 

end Fast_Magnetlc_Model; 
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B.l 4x4 Gravity Model Test Case from Ref. [2] 

Position vector = 5489150.0 . 802222.0 , 3140916.0 (meters) 

Gravity Model Data - Gem- 10 (See Ref (14J) 

_ 3 

4 = 398_600.47E9— - , r = 6_378_139.0 m. Degree and order = 4 
s 2 

•*** standard gravity model (As given In Ref. [21 )•*** Run Time: 2.70 Seconds 

Gravitational acceleration from acceleration only call 

-8.442692 12018857E+00 - 1 .23393633785485E+00 -4. 846593523466 14E+00 
Gravitational acceleration from acceleration plus dgdx 

-8.4426921 20 18857E+00 -1.23393633785485E+00 -4. 846593523466 14E+00 
Potential = 6.25359843440795E+07 
Analytic dgdx 

1.87779 19164782 IE-06 4.99270741320439E-07 1.96515882331 155E-06 
4. 9927074 1320439E-07 -1.46519995493325E-06 2.872141 12813307E-07 
1.96515882331 155E-06 2.872141 1281 3307E-07 -4.12591961544953E-07 

**** fast gravity model Run Time: 2.22 Seconds 

Gravitational acceleration from acceleration only call 

-8. 4426921 20 18857E+00 -1.23393633785485E+00 -4.84659352346614E+00 
Gravitational acceleration from acceleration plus dgdx 

-8. 442692 1201 8857E+00 -1.23393633785485E+00 -4.846593523466 14E+00 
Potential = 6.25359843440795E+07 
Analytic dgdx 

1 .87779 19164782 IE-06 4.9927074 1320439E-07 1.96515882331 155E-06 
4.9927074 1320439E-07 -1.46519995493325E-06 2.87214112813307E-07 
1.96515882331 155E-06 2.872141 12813307E-07 -4.12591961544953E-07 

**** normalized gravity model “norm_ir Run Time: 2.39 Seconds 

Gravitational acceleration from acceleration only call 

-8.442692 12018857E+00 -1.23393633785485E+00 -4. 846593523466 14E+00 
Gravitational acceleration from acceleration plus dgdx 

-8.442692 12018857E+00 -1.23393633785485E+00 -4.846593523466 14E+00 
Potential = 6.25359843440795E+07 
Analytic dgdx 

1 .87779 19164782 IE-06 4.99270741320439E-07 1.96515882331 155E-06 
4.9927074 1320439E-07 -1.46519995493325E-06 2.872 141 128 13307E-07 
1.96515882331 155E-06 2.872 141 128 13307E-07 -4.12591961544953E-07 

Note that the answers are identical to those in Ref. [2]. 
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B.2 5 z S Gravity Model Test Case from Ref. [2] 

Position vector = 5489150.0 . 802222.0 . 3140916.0 (meters) 

Gravity Model Data - Gem- 10 (See Ref [14)) 

3 

4 = 398 600.47E9-jr . r e = 6 378 139.0 m. Degree and order = 5 

s 2 

**** standard gravity model ( As given in Ref. [2] ) **** Run Time: 3.82 Seconds 

Gravitational acceleration from acceleration only call 

-8.44260633555472E+00 - 1.2339324305 1834E+00 -4.84652486332608E+00 
Gravitational acceleration from acceleration plus dgdx 

-8.44260633555472E+00 - 1.2339324305 1834E+00 -4.84652486332608E+00 
Potential = 6.25358693982450E+07 
Analytic dgdx 

1. 87773230503 190E-06 4.99259374934480E-07 1.965074721 12557E-06 
4.99259374934480E-07 -1.46513564895359E-06 2.8720884453 1796E-07 
1.965074721 12557E-06 2.8720884453 1796E-07 -4.12596656078305E-07 

•*** fast gravity model •••* Run Time: 3.17 Seconds 

Gravitational acceleration from acceleration only call 

-8.44260633555472E+00 - 1.2339324305 1834E+00 -4.84652486332608E+00 
Gravitational acceleration from acceleration plus dgdx 

-8.44260633555472E+00 -1. 2339324305 1834E+00 -4.84652486332608E+00 
Potential = 6.25358693982450E+07 
Analytic dgdx 

1. 87773230503 190E-06 4.99259374934480E-07 1.965074721 12557E-06 
4.99259374934480E-07 -1.46513564895359E-06 2.8720884453 1796E-07 
1.96507472 112557E-06 2.8720884453 1796E-07 -4.12596656078305E-07 

normalized gravity model “norm_H" **** Run Time: 3.38 Seconds 

Gravitational acceleration from acceleration only 

-8. 44260633555472E+00 -1. 2339324305 1834E+00 -4.84652486332608E+00 
Gravitational acceleration from acceleration plus dgdx 

-8. 44260633555472E+00 -1. 2339324305 1834E+00 -4.84652486332608E+00 
Potential = 6.25358693982450E+07 

1. 87773230503 190E-06 4.99259374934480E-07 1.96507472 112557E-06 
4.99259374934480E-07 -1.46513564895359E-06 2.8720884453 1796E-07 
1 .96507472 1 12557E-06 2.8720884453 1796E-07 -4. 12596656078305E-07 

Note that the answers are identical to those in Ref. [2). 
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B.3 Gravity Gradient Torque Teat Caae 

A comparison of general gravity gradient torque, computed with degree=2 order=0, 
with Rolthmyar ( Ref. 151 ) J2 only model is given below. 

Gravity Model - Gem- 10 (See Ref 1141) 

Position vector = 5489150.0 . 802222.0 .3140916.0 (meters) 

The rotation matrix from body to inertial is derived from 

pitch = 20.0 degrees 
yaw = 30.0 degrees 
roll = 40.0 degrees 

Mass tensor = ((477.0 , 63.0 , 0.0 ). 

(63.0 . 770.0, 0.0 ). 

(0.0 .0.0 .821.0)); 


gravity gradlent_torque from general gravity model (Spherical 0x0) 

-7.3839 160 15 19382E-05 -6.34664808264096E-04 3.51747050237606E-04 

gravity gradlent_torque from general gravity model (2x0) 
-7.35430183066330E-05 -6.34049979845033E-04 3.52179052151285E-04 

gravity gradlent_torque from Roithmayr J2 model 

-7.35430 18306633D-05 -6.3404997984503D-04 3.52 179052 15 128D-04 

gravity gradient_torque from general gravity model (4x4) 

-7.3549059240 1439E-05 -6.34 1023 12449005E-04 3. 522096 10370063E-04 

Note that the torque changes slightly as more terms are Included. 

The major impact Is expected to be control where the effect of 
forces that are longitude dependent will be included. 


58 


NASA CR- 188243 


January 1993 


Tir 


B.4 Magnetic Field Vector 

The position vector was (meters) 

X = (5489150.0 , 802222.0 , 3140916.0) 

Degree and Order = 10 

The resulting magnetic field vector was (Tesla) 

B = (-3.752 597530 1 8348E-05 , -6.1 600244200 1094E-06, 1.35121 1726546 19E-05) 
The mean radius of the earth was 6371.2 km 
The (10x10) harmonic data were taken from [1 1] Table II, IGRF 1985 
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Figure B.l 


Execution Tines for Fast, Normalized & Ref [2] Gravity Models 



This comparison show the relative efficiencies of the various implementations. Note that 
the normalized models, while less efficient than the unnormalized model, are none the 
less more efficient than the RefI2] model. This Is a result of the data structures and the 
precomputing of all possible derived Legendre functions. The simple normalized model 
(norm_I) and the normalized model (normjl), which uses the recursion relationship 
from [3J, have very similar speeds. In the next plot, norm_I Is not shown. Note also, that 
the “fast" model Is approximately 20% faster than the Ref[2] model. 
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Figure B.2 


Execution Tines for Fast, Normalized & Ref [2] Gravity Models 



This plot, which was generated using the 50x50 GEM-T3 model [151, shows relative effi- 
ciency as degree fit order increase. The difference in the run times compared to those in 
Fig B1 are a result of using a faster computer (Sun Sparc n Instead of Sun Sparc I). Note 
that past degree and order 27 the normalized model Is slower than the algorithm In [2J. 
This Is a result of the extra multiplication that must be done In the normalized methods 
finally overcoming the benefit of precomputing some of the data. The “fast" model Is 
always faster than [2]. It is worth noting that, at 50x50, and zero latitude, the gravita- 
tional acceleration vector computed using the simple normalized method (norm_I) dif- 
fered from the other three In the last three (out of 16) places. The benefit of the recursion 
relationship from (31 is beginning to make itself felt. The other three agreed to 15 places. 
The slight kinks in the plots are the result of timer noise and not some sudden change In 
computation speed. 
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